Offline Data Analysis of Strongyloides ratti bulk RNAseq data

Introduction —

The goal of this file is to pre-process the Strongyloides ratti RNAseq dataset originally published by Hunt et al 2016. Raw reads are aligned to the S. ratti reference transcription (PRJEB125.WBPS14.mRNA_transcripts, downloaded from WormBase Parasite on 17 August 2020)

Data Pre-Processing

A description of Kallisto alignment and data filtering/normalization steps can be found in Sr_RNAseq_Data_Processing.rmd.

Data Analysis

The limma package ( Ritchie et al 2015, Phipson et al 2016) is used to conduct pairwise differential gene expression analyses between life stages. Here, all unique pairwise comparisons between life stages are displayed as interactive volcano plots and DataTables. In this analysis, genes included in DataTables are corrected for multiple life-stage comparisons. An additional code chunk identifies all genes that are differentially expressed in any way across groups (i.e. an ANOVA analysis across life stages). This set of genes is passed into a clustering analysis that generates a heatmap of differential gene expression across life stages.

Load and Parse Preprocessed Data —

load (file = "../Outputs/SrRNAseq_data_preprocessed")
targets <- SrRNAseq.preprocessed.data$targets
annotations <- SrRNAseq.preprocessed.data$annotations
log2.cpm.filtered.norm <- SrRNAseq.preprocessed.data$log2.cpm.filtered.norm
myDGEList.filtered.norm <-SrRNAseq.preprocessed.data$myDGEList.filtered.norm

rm(SrRNAseq.preprocessed.data)

load(file = "../Strongyloides_RNAseq_Browser/Data/Sr_vDEGList")

Hierarchical Clustering and Principle Components Analysis —

# Introduction to this chunk -----------
# This code chunk starts with filtered and normalized abundance data in a data frame (not tidy).
# It will implement hierarchical clustering and PCA analyses on the data.
# It will plot various graphs and can save them in PDF files.
# Load packages ------
suppressPackageStartupMessages({
  library(tidyverse) # you're familiar with this fromt the past two lectures
  library(ggplot2)
  library(RColorBrewer)
  library(ggdendro)
  library(magrittr)
  library(factoextra)
  library(gridExtra)
  library(cowplot)
  library(dendextend)
})

# Identify variables of interest in study design file ----
group <- factor(targets$group)
batch <- factor(targets$block)
source <- factor(targets$source)

# Hierarchical clustering ---------------
# Remember: hierarchical clustering can only work on a data matrix, not a data frame

# Calculate distance matrix
# dist calculates distance between rows, so transpose data so that we get distance between samples.
# how similar are samples from each other
colnames(log2.cpm.filtered.norm)<-targets$group
distance <- dist(t(log2.cpm.filtered.norm), method = "maximum") #other distance methods are "euclidean", maximum", "manhattan", "canberra", "binary" or "minkowski"

# Calculate clusters to visualize differences. This is the hierarchical clustering.
# The methods here include: single (i.e. "friends-of-friends"), complete (i.e. complete linkage), and average (i.e. UPGMA). Here's a comparison of different types: https://en.wikipedia.org/wiki/UPGMA#Comparison_with_other_linkages
clusters <- hclust(distance, method = "complete") #other agglomeration methods are "ward.D", "ward.D2", "single", "complete", "average", "mcquitty", "median", or "centroid"
dend <- as.dendrogram(clusters) 

p1<-dend %>% 
  dendextend::set("branches_k_color", k = 5) %>% 
  dendextend::set("hang_leaves", c(0.1)) %>% 
  dendextend::set("labels_cex", c(0.5)) %>%
  dendextend::set("labels_colors", k = 5) %>% 
  dendextend::set("branches_lwd", c(0.7)) %>% 
  
  as.ggdend %>%
  ggplot (offset_labels = -0.5) +
  theme_dendro() +
  ylim(0, max(get_branches_heights(dend))) +
  labs(title = "Hierarchical Cluster Dendrogram ",
       subtitle = "filtered, TMM normalized",
       y = "Distance",
       x = "Life stage") +
  coord_fixed(1/2) +
  theme(axis.title.x = element_text(color = "black"),
        axis.title.y = element_text(angle = 90),
        axis.text.y = element_text(angle = 0),
        axis.line.y = element_line(color = "black"),
        axis.ticks.y = element_line(color = "black"),
        axis.ticks.length.y = unit(2, "mm"))
p1

# Principal component analysis (PCA) -------------
# this also works on a data matrix, not a data frame
pca.res <- prcomp(t(log2.cpm.filtered.norm), scale.=F, retx=T)
summary(pca.res) # Prints variance summary for all principal components.

Importance of components:
                            PC1     PC2     PC3      PC4      PC5     PC6     PC7     PC8
Standard deviation     113.3816 83.9229 76.8171 14.59647 12.62554 9.82530 8.57434 7.55633
Proportion of Variance   0.4833  0.2648  0.2219  0.00801  0.00599 0.00363 0.00276 0.00215
Cumulative Proportion    0.4833  0.7481  0.9700  0.97800  0.98400 0.98763 0.99039 0.99254
                           PC9    PC10    PC11    PC12    PC13      PC14
Standard deviation     7.23685 6.63777 6.15890 5.82959 5.48974 1.599e-13
Proportion of Variance 0.00197 0.00166 0.00143 0.00128 0.00113 0.000e+00
Cumulative Proportion  0.99451 0.99616 0.99759 0.99887 1.00000 1.000e+00

#pca.res$rotation #$rotation shows you how much each gene influenced each PC (called 'scores')
#pca.res$x # 'x' shows you how much each sample influenced each PC (called 'loadings')
#note that these have a magnitude and a direction (this is the basis for making a PCA plot)
## This generates a screeplot: a standard way to view eigenvalues for each PCA. Shows the proportion of variance accounted for by each PC. Plotting only the first 10 dimensions.
p2<-fviz_eig(pca.res,
             barcolor = brewer.pal(8,"Pastel2")[8],
             barfill = brewer.pal(8,"Pastel2")[8],
             linecolor = "black",
             main = "Scree plot: proportion of variance accounted for by each principal component",
             ggtheme = theme_bw()) 
p2

pc.var<-pca.res$sdev^2 # sdev^2 captures these eigenvalues from the PCA result
pc.per<-round(pc.var/sum(pc.var)*100, 1) # we can then use these eigenvalues to calculate the percentage variance explained by each PC

# Visualize the PCA result ------------------
#lets first plot any two PCs against each other
#We know how much each sample contributes to each PC (loadings), so let's plot
pca.res.df <- as_tibble(pca.res$x)

# Plotting PC1 and PC2
p3<-ggplot(pca.res.df) +
  aes(x=PC1, y=PC2, label=targets$group, 
      fill = targets$group,
      color = targets$group
  ) +
  geom_point(size=4, shape= 21, color = "black", alpha = 0.5) +
  #geom_label(color = "black", size = 2) +
  scale_fill_brewer(palette = "Set2") +
  scale_color_brewer(palette = "Set2", guide = FALSE) +
  #stat_ellipse() +
  xlab(paste0("PC1 (",pc.per[1],"%",")")) + 
  ylab(paste0("PC2 (",pc.per[2],"%",")")) +
  labs(title="Principal Components Analysis of S. stercoralis RNAseq Samples",
       sub = "Note: analysis is blind to life stage identity.",
       fill = "Life Stage") +
  scale_x_continuous(expand = c(.3, .3)) +
  scale_y_continuous(expand = c(.3, .3)) +
  coord_fixed() +
  theme_bw()
p3

# Create a PCA 'small multiples' chart ----
pca.res.df <- pca.res$x[,1:3] %>% 
  as_tibble() %>%
  add_column(sample = targets$sample,
             source = source,
             group = group,
             batch = factor(targets$block))

pca.pivot <- pivot_longer(pca.res.df, # dataframe to be pivoted
                          cols = PC1:PC3, # column names to be stored as a SINGLE variable
                          names_to = "PC", # name of that new variable (column)
                          values_to = "loadings") # name of new variable (column) storing all the values (data)
PC1<-subset(pca.pivot, PC == "PC1")
PC2 <-subset(pca.pivot, PC == "PC2")
#PC3 <- subset(pca.pivot, PC == "PC3")
#PC4 <- subset(pca.pivot, PC == "PC4")

p6<-ggplot(pca.pivot) +
  aes(x=sample, y=loadings) + # you could iteratively 'paint' different covariates onto this plot using the 'fill' aes
  geom_bar(stat="identity", fill = brewer.pal(8,"Pastel2")[8]) +
  scale_fill_brewer(palette = "Set2") +
  facet_wrap(~PC) +
  geom_bar(data = PC1, stat = "identity", aes(fill = batch)) +
  geom_bar(data = PC2, stat = "identity", aes(fill = batch)) +
  labs(title="PCA 'small multiples' plot",
       fill = "Life Stage Groups",
       caption=paste0("produced on ", Sys.time())) +
  scale_x_discrete(limits = targets$sample, labels = targets$group) +
  theme_bw() +
  coord_flip()
p6

# Graph all the plots generated in this script ----
# Plot all the plots on a grid
# p7<- grid.arrange(p1,p2,p3,p4, p5,p6,ncol = 4,
#             layout_matrix = cbind(c(1,1,1,2,2,2),c(3,3,4,4,5,5),c(6,6,6,6,6,6), c(6,6,6,6,6,6)))
# ggsave("Strongyloides stercoralis RNAseq Multivariate Analysis.pdf", 
#        plot = p7, 
#        device = "pdf",
#        #height = 17,
#        width = 22,
#        path = '../Outputs/'

Heatmap of Gene Expression Across Life Stages —

library(pheatmap)
library(RColorBrewer)
# Make a heatmap for all the genes using the Log2CPM values

diffGenes <- v.DEGList.filtered.norm$E %>%
  as_tibble(rownames = "geneID", .name_repair = "unique") %>%
  dplyr::select(!geneID) %>%
  as.matrix()
rownames(diffGenes) <- rownames(v.DEGList.filtered.norm$E)
colnames(diffGenes) <- as.character(v.DEGList.filtered.norm$targets$group)
clustColumns <- hclust(as.dist(1-cor(diffGenes, method="spearman")), method="complete")
clustRows <- hclust(as.dist(1-cor(t(diffGenes), 
                                  method="pearson")), 
                    method="complete") 
par(cex.main=1.2)

showticklabels <- c(TRUE,FALSE)
p<-pheatmap(diffGenes,
            color = RdBu(75),
            cluster_rows = clustRows,
            cluster_cols = clustColumns,
            show_rownames = F,
            scale = "row",
            angle_col = 45,
            main = "Log2 Counts Per Million (CPM) Expression Across Life Stages"
          
)

Differentially Expressed Genes —

# Introduction to this chunk ----
# This chunk uses a variance-stabilized DGEList of filtered and normalized abundance data.
# 
# These data/results are examples, a responsive version of this code is avaliable in a Shiny App.
# 
# Because we have access to biological and technical replicates, we can use statistical tools for differential expression analysis
# Useful reading on differential expression: https://ucdavis-bioinformatics-training.github.io/2018-June-RNA-Seq-Workshop/thursday/DE.html

# Load packages ----
suppressPackageStartupMessages({
  library(tidyverse)
  library(limma) # differential gene expression using linear modeling
  library(edgeR)
  library(gt) 
  library(DT) 
  library(plotly)
  library(ggthemes)
  library(RColorBrewer)
  source("../Strongyloides_RNAseq_Browser/Server/theme_Publication.R")
})

diffGenes.df <- v.DEGList.filtered.norm$E %>%
  as_tibble(rownames = "geneID", .name_repair = "unique")

# Set Expression threshold values for plotting and saving DEGs ----
adj.P.thresh <- 0.05
lfc.thresh <- 1 

group <- factor(v.DEGList.filtered.norm$targets$group)
design <- model.matrix(~0 + group) # no intercept/blocking for matrix, comparisons across group
colnames(design) <- levels(group)


# Fit a linear model to the data ----
fit <- lmFit(v.DEGList.filtered.norm, design = design)

# As an example, generate comparison matrix for a pairwise comparison ----
# iL3s vs FLF
# Note that the target/contrast goups will be divided by the number of life 
# stage groups e.g. PF+FLF/2 - iL3+iL3a+pfL1+ppL1+ppL3/5
comparison <- c('(iL3)-(FLF)')

targetStage<- comparison %>%
  str_split(pattern="-", simplify = T) %>%
  .[,1] %>%
  gsub("(", "", ., fixed = TRUE) %>%
  gsub(")", "", ., fixed = TRUE) %>%
  str_split(pattern = "\\+", simplify = T)

contrastStage<-comparison %>%
  str_split(pattern="-", simplify = T) %>%
  .[,2] %>%
  gsub("(", "", ., fixed = TRUE) %>%
  gsub(")", "", ., fixed = TRUE)  %>%
  str_split(pattern = "\\+", simplify = T)

comparison<- sapply(seq_along(comparison),function(x){
  tS <- as.vector(targetStage[x,]) %>%
    .[. != ""] 
  cS <- as.vector(contrastStage[x,]) %>%
    .[. != ""] 
  paste(paste0(tS, 
               collapse = "+") %>%
          paste0("(",.,")/",length(tS)),
        paste0(cS, 
               collapse = "+") %>%
          paste0("(",.,")/",length(cS)),
        sep = "-")
  
})

# Generate contrast matrix ----
contrast.matrix <- makeContrasts(contrasts = comparison,
                                 levels=design)

# extract the linear model fit -----
fits <- contrasts.fit(fit, contrast.matrix)
# empirical bayes smoothing of gene-wise standard deviations provides increased power (see: https://www.degruyter.com/doi/10.2202/1544-6115.1027)
ebFit <- eBayes(fits)

# Pull out the DEGs that pass a specific threshold for all pairwise comparisons ----
# Adjust for multiple comparisons using method = global. 
results <- decideTests(ebFit, method="global", adjust.method="BH", p.value = adj.P.thresh)

recode01<- function(x){
  case_when(x == 1 ~ "Up",
            x == -1 ~ "Down",
            x == 0 ~ "NotSig")
}
diffDesc <- results %>%
  as_tibble(rownames = "geneID") %>%
  dplyr::mutate(across(-geneID, unclass)) %>%
  dplyr::mutate(across(where(is.double), recode01))

# Function that identifies top DEGs between a specific contrast ----
calc_DEG_tbl <- function (ebFit, coef) {
  myTopHits.df <- limma::topTable(ebFit, adjust ="BH", 
                                  coef=coef, number=40000, 
                                  sort.by="logFC") %>%
    as_tibble(rownames = "geneID") %>%
    dplyr::rename(tStatistic = t, LogOdds = B, BH.adj.P.Val = adj.P.Val) %>%
    dplyr::relocate(UniProtKB, Description, InterPro, GO_term, WBgeneID, Str_geneID, Str_percent_homology, Ce_geneID, Ce_percent_homology, .after = LogOdds)
  
  myTopHits.df
}

list.myTopHits.df <- sapply(comparison, function(y){
  calc_DEG_tbl(ebFit, y)}, 
  simplify = FALSE, 
  USE.NAMES = TRUE)

list.myTopHits.df <- sapply(comparison, function(y){
  list.myTopHits.df[[y]] %>%
    dplyr::select(geneID, 
                  logFC, 
                  BH.adj.P.Val:Ce_percent_homology)},
  simplify = FALSE, 
  USE.NAMES = TRUE)

# Get log2CPM values and threshold information for genes of interest
list.myTopHits.df <- sapply(seq_along(comparison), function(y){
  tS<- targetStage[y,][targetStage[y,]!=""]
  cS<- contrastStage[y,][contrastStage[y,]!=""]
  
  concat_name <- function(x) {
    ifelse(x == "target", 
           paste(tS, collapse = "+"), 
           paste(cS, collapse = "+"))
  }
  
  groupAvgs <- diffGenes.df %>%
    dplyr::select(geneID, starts_with(paste0(tS,"-")), 
                  starts_with(paste0(cS,"-"))) %>%
    pivot_longer(cols = -geneID, names_to = c("group","sample"), values_to = "CPM",
                 names_sep = "-") %>%
    dplyr::mutate(contrastID = if_else(group %in% tS,"target", "contrast")) %>%
    group_by(geneID, contrastID) %>%
    dplyr::select(-sample) %>%
    summarize(mean = mean(CPM), .groups = "drop_last") %>%
    pivot_wider(names_from = contrastID, values_from = mean) %>%
    dplyr::relocate(contrast, .after = target) %>%
    dplyr::rename_with(concat_name, -geneID) %>%
    dplyr::rename_with(.cols =-geneID, .fn = ~ paste0("avg_(",.x,")"))
  
  diffGenes.df %>%
    dplyr::select(geneID, starts_with(paste0(tS,"-")), 
                  starts_with(paste0(cS,"-"))) %>%
    left_join(groupAvgs, by = "geneID") %>%
    left_join(list.myTopHits.df[[y]],., by = "geneID") %>%
    left_join(dplyr::select(diffDesc,geneID,comparison[y]), by = "geneID") %>%
    dplyr::rename(DEG_Desc=comparison[y]) %>%
    dplyr::relocate(DEG_Desc) %>%
    dplyr::relocate(logFC:Ce_percent_homology, .after = last_col())
  
},
simplify = FALSE)

comparison <- gsub("/[0-9]*","", comparison)
names(list.myTopHits.df) <- comparison

list.myTopHits.df <- sapply(comparison, function(y){
  list.myTopHits.df[[y]] %>%
    dplyr::mutate(DEG_Desc = case_when(DEG_Desc == "Up" ~ paste0("Up in ", str_split(y,'-',simplify = T)[1,1]),
                                       DEG_Desc == "Down" ~ paste0("Down in ", str_split(y,'-',simplify = T)[1,1]),
                                       DEG_Desc == "NotSig" ~ "NotSig")) 
},
simplify = FALSE, 
USE.NAMES = TRUE)

# PC1 Volcano Plot and Interactive Table ----
vplot1 <- ggplot(list.myTopHits.df[[1]]) +
  aes(y=-log10(BH.adj.P.Val), x=logFC, text = paste(geneID, "<br>",
                                                    "logFC:", round(logFC, digits = 2), "<br>",
                                                    "p-val:", format(BH.adj.P.Val, digits = 3, scientific = TRUE))) +
  geom_point(size=2) +
  geom_hline(yintercept = -log10(adj.P.thresh), 
             linetype="longdash", 
             colour="grey", 
             size=1) + 
  geom_vline(xintercept = lfc.thresh, 
             linetype="longdash", 
             colour="#BE684D", 
             size=1) +
  geom_vline(xintercept = -lfc.thresh, 
             linetype="longdash", 
             colour="#2C467A", 
             size=1) +
  labs(title = paste0('Pairwise Comparison: ',
                      gsub('-',
                           ' vs ',
                           comparison[1])),
       subtitle = paste0("grey line: p = ",
                         adj.P.thresh, "; colored lines: log-fold change = ", lfc.thresh),
       color = "GeneIDs") +
  theme_Publication() 
vplot1

# Interactive Tables
yy<- 1
tS<- targetStage[yy,][targetStage[yy,]!=""]
cS<- contrastStage[yy,][contrastStage[yy,]!=""]
sample.num.tS <- sapply(tS, function(x) {colSums(v.DEGList.filtered.norm$design)[[x]]}) %>% sum()
sample.num.cS <- sapply(cS, function(x) {colSums(v.DEGList.filtered.norm$design)[[x]]}) %>% sum()


n_num_cols <- sample.num.tS + sample.num.cS + 5


LS.datatable <- list.myTopHits.df[[yy]] %>%
  DT::datatable(rownames = FALSE,
                caption = htmltools::tags$caption(
                  style = 'caption-side: top; text-align: left; color: black',
                  htmltools::tags$b('Differentially Expressed Genes in', 
                                    htmltools::tags$em('S. ratti'), 
                                    gsub('-',' vs ',comparison[yy])),
                  htmltools::tags$br(),
                  "Threshold: p < ",
                  adj.P.thresh, "; log-fold change > ",
                  lfc.thresh,
                  htmltools::tags$br(),
                  'Values = log2 counts per million'),
                options = list(autoWidth = TRUE,
                               scrollX = TRUE,
                               scrollY = '300px',
                               scrollCollapse = TRUE,
                               order = list(n_num_cols-1, 
                                            'desc'),
                               searchHighlight = TRUE, 
                               pageLength = 25, 
                               lengthMenu = c("5",
                                              "10",
                                              "25",
                                              "50",
                                              "100"),
                               columnDefs = list(
                                 list(
                                   targets = ((n_num_cols + 
                                                 1)),
                                   render = JS(
                                     "function(data, row) {",
                                     "data.toExponential(1);",
                                     "}")
                                 ),
                                 list(
                                   targets = ((n_num_cols + 
                                                 4):(n_num_cols + 
                                                       5)),
                                   render = JS(
                                     "function(data, type, row, meta) {",
                                     "return type === 'display' && data.length > 20 ?",
                                     "'<span title=\"' + data + '\">' + data.substr(0, 20) + '...</span>' : data;",
                                     "}")
                                 ),
                                 list(targets = "_all",
                                      class="dt-right")
                               ),
                               rowCallback = JS(c(
                                 "function(row, data){",
                                 "  for(var i=0; i<data.length; i++){",
                                 "    if(data[i] === null){",
                                 "      $('td:eq('+i+')', row).html('NA')",
                                 "        .css({'color': 'rgb(151,151,151)', 'font-style': 'italic'});",
                                 "    }",
                                 "  }",
                                 "}"  
                               ))
                               
                )) 
LS.datatable <- LS.datatable %>%
  DT::formatRound(columns=c(3:n_num_cols), 
                  digits=3)

LS.datatable <- LS.datatable %>%
  DT::formatRound(columns=c(n_num_cols+2, 
                            n_num_cols+9,
                            n_num_cols+11), 
                  digits=2)

LS.datatable <- LS.datatable %>%
  DT::formatSignif(columns=c(n_num_cols+1), 
                   digits=3)

LS.datatable

NA

Benchmarking

suppressPackageStartupMessages({
  library(openxlsx)
  library(tidyverse)
  library(ggplot2)
})
# Load Hunt Dataset: iL3 vs FLF comparison
temp.dat <-  read.xlsx ("../Benchmarking/ng.3495-S5.xlsx", 
                        sheet = 2, startRow = 2)

Hunt.dat <- tibble(geneID = temp.dat$X2, logFC = temp.dat$logFC)
rm(temp.dat)

# Rename Results of iL3 vs FLF comparison from Browser
Browser.dat <- list.myTopHits.df$`(iL3)-(FLF)` %>%
  dplyr::select(geneID, logFC)

print(paste('Total number of genes in Hunt *et al* 2016 iL3 vs FLF comparison tab:',nrow(Hunt.dat)))

[1] "Total number of genes in Hunt *et al* 2016 iL3 vs FLF comparison tab: 11146"

print(paste('Total number of genes in Str-RNAseq Browser iL3 vs FLF output file:', nrow(Browser.dat))) 

[1] "Total number of genes in Str-RNAseq Browser iL3 vs FLF output file: 12035"

# The plot below takes the genes with LogFC results in both the Browser and Hunt databases, and plots the two sets against each other. 
plotting.all <- inner_join(Browser.dat, Hunt.dat, by = "geneID")
ggplot(plotting.all, aes(x = logFC.x, y = logFC.y)) +
  geom_smooth(method=lm, color = 'red', formula = "y ~ x") +
  geom_point(shape=16, size=3, alpha = 0.8) +
  labs(title = "All Overlapping Genes: Str-Browser vs Hunt Data",
       subtitle = "iL3 vs FLF",
       caption = "points = genes; red line = linear model (y ~ x)",
       x = "Str-Browser LogFC",
       y = "Hunt LogFC") +
  coord_equal() +
  theme_bw()

NA
NA

Cluster DEGs into functional modules —

Functional Enrichment Analysis —

nrow(GO.present.Compara.absent)

[1] 543